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Abstract 

In this paper the global eigenmode structures of linear ion-temperature-gradient (ITG) modes in 
tokamak plasmas are obtained using a novel technique which combines results from the local gyrokinetic 
code GS 2 with analytical theory to reconstruct global properties. Local gyrokinetic calculations, using 
GS 2 , are performed for a range of radial flux surfaces, x , and ballooning phase angles, p , to map out the 
local complex mode frequency, Qq(x,p) = uJo(x,p)+i'yo(x,p) for a single toroidal mode number, n. Taylor 
expanding about a reference surface at x = 0, and employing the Fourier-ballooning representation 
leads to a second order ODE for the amplitude envelope, A (p), which describes how the local results are 
combined to form the global mode. The equilibrium profiles impact on the variation of Oo (x,p) and hence 
influence the global mode structure. The simulations presented here are based upon a global extension 
to the CYCLONE base case and employ the circular Miller equilibrium model. In an equilibrium with 
radially varying profiles of a/L^ and a/L n , peaked at x = 0, and with all other equilibrium profiles held 
constant, including r]i = L u /Lt , £lo(x,p) is found to have a stationary point. The reconstructed global 
mode sits at the outboard mid-plane of the tokamak, with global growth rate, 7 ~Max[ 7 0 ]. Including the 
radial variation of other equilibrium profiles like safety factor and magnetic shear, leads to a mode that 
peaks away from the outboard mid-plane, with a reduced global growth rate. Finally, the influence of 
toroidal flow shear has also been investigated through the introduction of a Doppler shift, cco —> uJo—ntl'^x, 
where W is the equilibrium toroidal flow, and a prime denotes the radial derivative. The equilibrium 
profile variations introduce an asymmetry into the global growth rate spectrum with respect to the sign 
of such that the maximum growth rate is achieved with non-zero shearing, consistent with recent 
global gyrokinetic calculations. 


1 Introduction 

Tokamaks [1] provide one of the most stable and promising configurations for magnetic confinement fusion. 
However, confinement in tokamaks is not perfect; there are a number of mechanisms by which energy 
and particles can be transported across the magnetic flux surfaces from the core confinement region to 
the plasma edge. The main contribution is typically due to turbulent transport, which is widely believed 
to originate primarily from microinstabilities driven by density and temperature gradients. These “drift” 
modes, with low frequency compared to the cyclotron frequency, are the dominant tokamak microinstabilities 
[1,2], and the turbulence they drive influences the minimum size of magnetic confinement fusion reactors, 
such as ITER and DEMO [3]. It is therefore important to understand these drift modes in order to seek 
ways to reduce their impact. Previous theoretical and numerical studies have shown that flow shear can 
control the stability of drift waves, providing a mechanism to suppress or even stabilise them completely 
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[4-9]. While microinstabilities must ultimately be saturated by nonlinear physics (e.g. zonal flows), it 
nevertheless remains important to understand the structure of linear instabilities and the threshold gradients 
associated with their onset: linear physics lies close to the heart of many simplified plasma models that are 
of considerable value, e.g. the quasi-linear based TGLF model of anomalous transport in tokamaks [10]. 
Microinstabilities are often investigated via numerical solutions of the gyrokinetic equation [11-13] and 
several different approaches are typically used, as we now discuss. 

For high toroidal mode numbers, n, the rational flux surfaces (i.e. those where the safety factor, q , is 
rational) are closely packed so equilibrium quantities vary only weakly from one surface to the next. Two 
length scales can then be identified: the equilibrium length scale, characterised by the minor radius, a, and 
the distance between rational surfaces, A = (n dq/dr ) _1 , where the derivative is with respect to radius, r. 
The separation between these length scales is exploited in ballooning theory [14-17] to perform an expansion 
in the small parameter, A/a. The rational surfaces spanned by a mode are then equivalent to leading order. 
Local gyrokinetic codes, like GS2 [16,18], exploit this “ballooning symmetry” to reduce the intrinsic 2D 
problem (in radius, r, and poloidal angle, 6) to ID in the extended ballooning coordinate, 77 , which is 
aligned with magnetic field lines. Local gyrokinetic codes only provide the structure along the magnetic 
field line, together with the local eigenvalue, flo(x,p) (where x = (r — ro) /a is the normalised radial distance 
from the reference surface, ro, and the ballooning phase angle, p, is the value of 77 where the radial derivative 
of the eigenfunction is zero). Local codes, however, provide neither the radial mode structure nor the global 
eigenvalue Q, as both x and p are free parameters at this lowest order in the ballooning expansion where 
local gyrokinetics is valid. The radial mode structure and global eigenvalue, D, are determined at the next 
order in the A/a expansion, where the eigenfunction dependence on x and p becomes constrained. 

Exploiting the higher order theory to solve for the full global eignmode has been demonstrated previously 
for a simplified fluid model of ion temperature gradient (ITG) modes [19]. In this paper, we build on Ref [19] 
to show how this formalism can be extended to more realistic full gyrokinetic plasma models. Our approach 
exploits the GS 2 code [16,18] to provide the local mode structure, £(z,p, 77 ), and Do(x,p), from which the 
higher order theory provides the global mode structure and D. 

It has previously been demonstrated that under the very special conditions where flo(x,p) has a station¬ 
ary point, a so-called isolated or conventional ballooning mode is obtained from the higher order analysis [14]. 
This type of mode, originally studied in the context of ideal MHD theory [20,21], is usually located at the 
outboard mid-plane, where the poloidal angle 9 = 0 [19]. The higher order theory for such an isolated mode 
provides a global complex mode frequency D, which includes an 0(l/n) correction to the local eigenvalue, 
Do(x,p), evaluated at the maximally unstable x and p. This maximally unstable ”isolated mode” requires 
a stationary point in the local complex mode frequency Qo(x,p), which is a very special situation. More 
usually a second class of mode, known as the “general mode”, would be expected in most situations. The 
higher order theory for these general modes predicts that they peak away from the outboard mid-plane and 
have a reduced growth rate relative to the isolated mode. Local gyrokinetic simulations alone cannot distin- 
gush between these two types of mode, but evidence for both can be found in global gyrokinetic simulations 
that include the radial variation of equilibrium profiles. For example, electrostatic ITG modes are found to 
be shifted relative to the outboard mid-plane in linear global gyrokinetic simulations of ASDEX Upgrade 
plasmas [ 22 ]. Such up-down asymmetries are important as they could provide a mechanism for flow gener¬ 
ation in tokamaks [23-26], which may be important in a machine like ITER for which the external torque 
is small. While nonlinear simulations are likely neccessary for a complete understanding of turbulence and 
flows in plasmas, linear theory provides a picture of the important physical mechanisms. The technique we 
present here is an alternative approach to full global simulations that uses an efficient formalism to shed 
light on the key linear physics. 

The paper is organised as follows. In Section 2 we review the theoretical formalism on which this paper 
is based. Results from applying this technique to the so-called CYCLONE base case [27,28], are presented 
in Section 3. Employing this standard test case enables us to compare with previously published global 
simulations [29]. Finally, Section 4 summarises our conclusions, and plans for future work. 


2 


2 The technique: From local to global gyrokinetic calculations 


We employ the initial value local gyrokinetic code, GS 2 [16,18], to solve the linearised gyrokinetic equation 
[ 11 , 12 ] numerically. GS 2 provides the local eigenvalue, flo(x,p) as well as the mode structure along a given 
magnetic field line £(a;,p, 77) in the infinite domain in 77 with a periodic dependence on p. To reconstruct the 
linear global mode properties from these local modes, we employ the Fourier-Ballooning (FB) representation 
[30]: 

/ OO 

£(:r,p, #, t) exp(— inq^O) exp \—inq'x(6 — p)\A(p)dp ( 1 ) 

-00 

where (j)(x,6,t ) is the global mode structure for fluctuations in the electrostatic potential, x is the radial 
variable, go is the safety factor at x = 0, and the local mode structure £(:r,p, #,£) is obtained from GS 2 . 
Note that, the ^-dependence in £ is due to a slow dependence of the equilibrium on x , and is a parameter 
at this order. For the linear global modes studied in this paper, £>(a;,#,t) and £(x,p, 6 A) have a separable 
time dependence of the form oc e~ lQlt and e ~ z ^o(x,p)t w here ft and flo(x,p) are the global and local complex 
linear eigenmode frequencies, respectively. In the rest of the paper, we will frequently drop the explicit 
t dependence, and use <j)(x,6) and £(x,p,0) to denote the separated spatial dependent part of the linear 
eigenmode structure. The mapping in Eq. (1) from the infinite domain in 77 to the periodic poloidal angle 
6 , is possible because of the symmetry property £(x,p + 2/7T, 6 + 2hr) = £(x,p, 0 ), for any integer l. 

In order to employ Eq. (1) to find the global mode structure, 0(a;, 0,t), we must evaluate the envelope 
A(p), which can be obtained from the higher order theory as follows. We seek a global eigenmode satisfying 


d<j>(x, 6 , t) 
dt 


—iQ</)(x, 6,t) 


( 2 ) 


where ft = uj + 77 , with uj and 7 corresponding to the global frequency and growth rate, respectively. 
Substituting Eq. (1) into Eq. (2), and noting d£(x,p,9,t)/dt = —ifto(x,p)^(x,p,6,t), gives 


/ OO 

-00 


fto(x,p)f;(x,p, 0, t) exp(— inqo6) exp [—inq'x(0 — p)\A(p)dp 


(3) 


Using Eq. (1) for </>, we can then write our eigenmode condition in the form: 



[ft — flo(x,p)] £(a;,p, 0 , t) ex.p(—inqo6) exp [—inq'x(6 — p)\A(p)dp = 0 


(4) 


Here, the local complex mode frequency, fto(x,p) and £(z,p, 0, £) are mapped out by running the local 
gyrokinetic code, GS 2 , several times spanning the required range of x and over — n < p < n. This process is 
trivially parallelised, and therefore will usually be more efficient than a full global simulation. Anticipating 
radially localised solutions in x we may Taylor expand fto to second order in x about x = 0 and Fourier 
expand in p to write: 

N 2 

n 0 (x,p) = ^2Y1 4 m) x m cos{kp) (5) 

k=0 m =0 

where N is the number of Fourier modes retained. The coefficients a are complex numbers, obtained by 
fitting the expansion to the full fto(x,p) results which are obtained from GS 2 . Now substituting Eq. (5) 
into Eq. (4), we obtain 



n - SkLo (4 0) +4U+4V) cos(fcp)] A(p) 

£(a;,p, 6A) exp(— inqoO) exp [—inq f x(6 — p)\dp = 0 . 


We note that nq'p is a radial wave number. Therefore, following standard Fourier transform procedures and 
assuming A(p) varies much more rapidly than £ withp, allows us to replace x m A(p) with (—i/nq') rn d m A/dp rn , 
to give: 



&A - ££Lo ( cos (kp) - Tjr cos (kp)± - cos (kp)^ ) A 


,(o) 


nq' 


d_ 

dp 


4 2) 

(nq') 2 


d 2 


£(a;,p, 6A) exp(— inqoO) exp [—inq'x(6 — p)\dp = 0. 


( 6 ) 
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Figure 1: From left to right: The variation of (a) real frequency, c^o, (b) linear growth rate, 70 , with k y pi, 
for the dominant modes at p = 0 for two values of the normalised ion collision frequency, vna/c s = 0 (solid 
line) and vna/c s = 0.28 (solid line with * symbols); (c) shows the local mode structure, £(x = 0,p,p = 0 ), 
for the most unstable mode at k y pi = 0.58, as a function of ballooning coordinate 77 along the magnetic field 
line. Note that both (jJq and 70 are measured in units of ( c s /a ) and these local simulations have been carried 
out at mid-radius, i.e. r = a/2. 


Equation ( 6 ) must hold for all x and 0 , which then provides our final equation for A(jp): 


N k 

E 

k =0 L 


,( 2 ) 


;J1) 


ia 


C 0 S (M-E - 7 E cos (W t + 4 0) cos {kp) 


(nq') 2 


dp 


nq 1 


dp 


A = ttA 


(7) 


As <fi(x,0,t) must be periodic in 6 , the representation Eq. ( 1 ) requires A{p) to be periodic in p. Thus Eq. (7) 
must be solved with periodic boundary conditions to determine the envelope A(p) and the global complex 
mode frequency, D as an eigenvalue. Knowledge of A(p), together with £(x,p,0) from GS 2 then allows the 
full 2D eigenfunction, </>(#, 0), to be reconstructed from Eq. (1). 

Equation (7) can be solved analytically in the limiting cases where either = 0, or where = 0 for 
all k ; the former implies that Qq(x,p) is stationary at x = 0, i.e. dQj$/dx\ x= o = 0. These two limits lead 
to two different classes of eigenmode [14,19] which are referred to as “isolated modes” in the special case 
where = 0 , and “general modes” in the more usual situation where 7 ^ 0 (when the terms in afjp can 
be neglected in the limit of large n). In this paper we shall retain both terms and solve Eq. (7) numerically. 


3 Results 

Before considering the global calculations, we first describe the local ballooning analysis. Our focus is on 
ITG modes, believed to be one of the common causes for turbulent transport in many tokamak plasmas. 
Nevertheless, we note that our technique applies to all classes of high n toroidal drift modes. We adopt the 
CYCLONE parameters [27,28] and for convenience we have limited ourselves to electrostatic calculations 
with adiabatic electrons, which allows us to make comparisons with previously published results from global 
simulations. The model parameters are given in Table 1 for a Miller equilibrium with circular flux surfaces. 

Figure 1 shows the local real frequency, cjo, and growth rate, 70 , obtained from GS 2 as functions of 
normalised binormal wavenumber, k y pi, for the dominant modes with p = 0 , together with the local mode 
structure, £ at x = 0. The results for a normalised ion-ion collision frequency of vaa/c s — 0.28 are found to 
be very similar to those for a collisionless plasma for p = 0. We use this finite collision rate in the remaining 
calculations as it helps damp unphysical modes found by GS 2 at values of p close to marginal stability, 
and yet gives local eigenmode result similar to the collisionless case. The most unstable mode is found at 
k y pi = 0.58 which corresponds to toroidal mode number n = 39. Now we apply the technique presented in 
Section 2 to reconstruct the global structure for the most unstable mode. 

3.1 Global Calculations: Isolated modes with flat 77 * profile 

We first consider radially varying profiles of temperature and density, such that the gradients a/L^ and a/L n 
peak in the centre of the domain at x = 0 while maintaining a flat profile in vp = L u /Lt (see Figure 2c). 
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Figure 2: The radial profiles used in this pa¬ 
per, chosen to match Ref [29]. (a) Safety fac¬ 

tor, q = 0.84 + 2.24(r/a) 2 and magnetic shear 
s = 2 ^1 — and (b) temperature, 
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Table 1: CYCLONE equilibrium parameters. 


Parameter 

Value 

Parameter 

Value 

s 

0.8 

r 0 (m) 

0.313 

Qo 

1.4 

p 

0.0 

a / Lj' 

2.54 

nq' 

144 

a/L n 

0.81 

vua.lc s 

0.28 

ky Pi 

0.58 

T P 

1.0 

a(m) 

0.625 

Pi(m) 

0.003384 

R(m) 

1.70 

n — 

P*~ a 

0.005415 


Table 2: The model coefficients, a™, with ten Fourier modes. The real and imaginary components contribute 
to the real frequency, cjo, and linear growth rate, 70 , respectively. Note that coefficients with m = 1 are all 
zero for this special case where only profiles in a/L^ and a/L n are retained. 


m 

k 

0 

2 

0 

0.1177 - 0.0680 i 

-1.5689 - 1.9352 i 

1 

0.1804 + 0.1221 i 

1.3347 - 2.2466 i 

2 

0.0462 + 0.0461 i 

0.0825 - 0.8734 i 

3 

0.0229 + 0.0231 i 

-0.0015 - 0.2477 i 

4 

0.0068 + 0.0098 i 

-0.1007 - 0.0227 i 

5 

0.0012 + 0.0078 i 

-0.1090 + 0.1078 i 

6 

-0.0022 + 0.0045 i 

-0.1134 + 0.1240 i 

7 

-0.0033 + 0.0023 i 

-0.0861 + 0.0856 i 

8 

-0.0035 - 0.0000 i 

-0.0461 + 0.0105 i 

9 

-0.0026 - 0.0018 i 

0.0111 - 0.0525 i 


Initially we fix all other equilibrium profiles to the constant values given in Table 1 (ie. we use the dashed 
line profiles of Figure 2a and b) 1 . GS 2 calculations across the radial extent of the mode in x and full range 
in ballooning phase angle — n < p < n then provide the local complex mode frequency, Qq(x,p). Figure 3 
shows Qq(x,p) from GS2, compared to the fit using the parameterisation given in Eq. (5). The coefficients 
resulting from the fit are given in Table 2 . As expected, because the linear drive profiles peak and are 
symmetric about x = 0 , the coefficients with m = 1 are negligible and Uo has a stationary point at x = 0 , 

p = 0 . 

Solving Eq. (7) numerically with the values of a™ in Table 2 we obtain A(p) and the associated eigenvalue 
Q = 0.362 + 0.135b Using the numerical solution for A(p), and £(a;,p, 6) obtained from GS2, Eq. (1) provides 
the global mode structure, <j)(x^0). Figure 4 shows our solution for A(p) and the corresponding solution for 
6) in the poloidal cross-section. We see that A(p) is localised, as is required for the procedure to be 
accurate (recall, we assumed A(jp) varies rapidly with p to derive Eq. ( 6 )). As shown in Figure 4, this leads 
to a mode that balloons on the outboard mid-plane at 6 = 0 . 

The radial mode width, A w r , and its variation with the toroidal mode number, n (or equivalently 1/p*), 
has also been calculated and is shown in Figure 5. Here the width is defined as the full width half maximum 
of a Gaussian fit to the magnitude of 0 at 6 = po , where po is the ballooning phase angle at which the global 
mode peaks poloidally. We find that the radial width of the mode scales inversely with the square root of 
n, Aw r oc n _CU9 as expected for isolated modes [20]. 

Changing the toroidal mode number (and p* to keep k y pi fixed) also affects the global mode frequency, 
U. Figure 6 shows that both the real frequency, cu, and the linear growth rate, 7 , scale with n according to: 
oj = 0.364 — 0.072n~ 1-12 and 7 = 0.147 — 0.492n -1-02 , respectively. This indicates that the finite n correction 
scales inversely with the toroidal mode number as expected from conventional ballooning theory [20]. In 
the limit n —»• 00 , the global complex mode frequency U converges to the local complex mode frequency at 
x = p = 0, Uo(0, 0) = £Jo(0, 0) + 270 ( 0 , 0) = 0.364 + 0.147b This confirms that our results are consistent with 

1 Results including the equilibrium variation shown in Figure 2a and b will be considered in the next subsection. 
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Figure 3: Contour plots of real and imaginary parts of the local complex mode frequency, measured in unit 
of (c s /a), as functions of radius, x , and ballooning phase angle, p, for parabolic Lt and L n radial profiles 
(see Figure 2c), while excluding other profile variations, (a) and (b) are, respectively, the frequency and 
growth rate obtained from the local gyrokinetic code, GS 2 . The corresponding frequency and growth rate 
from the fitted model, using Eq. (5), are presented in (c) and (d) respectively. The * symbols indicate the 
marginal stability contour where 70 (x : p) = 0 . 



Figure 4: (a) presents the numerical solution of the envelope function, 4(p), obtained from Eq. (7), as a 
function of ballooning phase angle, p, and (b) shows the reconstructed electrostatic potential global mode 
structure, (j)(x , 0), in the poloidal plane. The profile variations other than Lt and L n are excluded from this 
calculation, and the mode structure, centred on the outboard mid-plane, is aligned radially where 6 = 0. 
The solid black lines indicate the radial domain of the calculation. 
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Figure 5: The radial mode width, Aw r , as a function 
of toroidal mode number, n. It scales, approximately, 
inversely with square root of n according to: A w r = 
0.687n -0 - 49 , as expected for an isolated mode. 


a b 




Figure 6 : (a) The global real frequency lj and (b) growth rate 7 as functions of toroidal mode number, n, for 
fixed kypi = 0.58. The uj and 7 are measured in units of ( c s /a ). They scale with n as: uj = 0.364—0.072n -1 ’ 12 , 
7 sb 0.147 — 0.492n -1 ’ 02 . 


analytic theory and higher order ballooning calculations presented in references [14,19,31]. 

To sumarise, this mode has all the characterstics of the isolated mode identified in [14,15,19]. It exists 
in this particular case because our choice of profiles ensure that f 2 o(^,p) has a stationary point at x == 0 , 
p = 0, as is clear from Figure 3. However, as we shall see in the following subsection, taking into account the 
radial variation of other equilibrium profiles introduces a small, but significant, deviation from conventional 
ballooning (or isolated) modes. 


3.2 Global Calculations with Profile Variations 

Now we introduce the other radially varying profiles given in Figure 2 , and repeat the analysis of subsection 
3.1. These profiles vary over the radial scale of the instability and give rise to a significant linear radial 
variation in Qq (x,p), such that the terms in Eq. (5) become significant. These extra linear terms impact 
ojo (x,p) and 70 (x,p) differently, which influences the global mode structure. Specifically, the positions of 
the maxima in ujq(x,p) and 70 (x, p) occur at different x , which affects the stability and poloidal position of 
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Figure 7: The reconstructed electrostatic potential global mode structure, 0), for n = 39, in the poloidal 
plane for different radial profile variations taken from Figure 2: (a) Lt, L n and x vary, here both temperature 
T and density n e are assumed to be constant, (b) Ly, L n , q and s vary (c) Ly, T n , x, q and s vary and 
finally (d) full profile variation in which Ly, T, L n , n e , x, q and s all vary. 




Figure 8 : The linear growth rate, 7 , measured in unit of (c s /a), as a function of flow shear, calculated 
for the most unstable mode with k y pi = 0.58. The toroidal mode number, n = 39 and q' « 3.6. For curve 
(e) the radial variation of the profiles other than Lt and L n are excluded, while the other curves correspond 
to the profile variations of Figure 7(a-d) respectively. 
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Figure 9: The scaling with toroidal mode numbeer, n, of (a) the value of flow shear to get the most unstable 
mode (crosses), fitted to the curve |f^ m | ~ 0.713n~ 1-03 , and (b) the offset in poloidal angle (crosses), fitted 
to the curve |po| — 0.041 + 1.88n -0-925 , for the full profile variation case, i.e curve—d from Figure 8. 

/ / / / 
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Figure 10: Reconstructed global mode structure, 0, in a small region of the poloidal cross-section at the 
outboard mid-plane (see Figure 4), under the combined effects of both radially varying equilibrium profiles 
from Figure 7-d and flow shear, From left to right, = —0.03,—0.017, 0 and 0.008, respectively. For 
= 0, the mode is already tilted, due to the profile variation effect, while a critical value of flow shear, 
occurs at « —0.017, which cancels out the effect of profile variation and, once again, the mode structure 
is aligned radially (as for a conventional ballooning mode). From Figure 8, we see that this is the maximally 
unstable flow shear. 


the global mode. 

Figure 7 shows the reconstructed global mode structure for a number of cases where we introduce 
different profiles. The effect of variation in radius x itself is presented in panel (a). The mode shifts slightly 
downward with respect to the outboard mid-plane. Panel (b) shows the effect of q and s profile variations 
and results in a mode with an upwards poloidal shift. Panel (c) shows the combined influence of both (a) 
and (b) and demonstrates that for this particular case they are competing and almost cancel. The result 
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is slightly in favour of the x variation, leading to a slight net downward poloidal shift with respect to the 
outboard mid-plane. Finally, in panel (d), in addition to the profiles considered in (c), we include radial 
variations in the T and n e profiles. The reconstructed global mode is then shifted poloidally downward with 
respect to the outboard mid-plane. 

For the various profile variations considered above, the up-down poloidal symmetry is broken. This 
moves the mode slightly off the mid-plane and towards the good curvature region, which reduces the growth 
rate compared to the flat r]i profile case considered in subsection 3.1. Specifically, for case d, with full profile 
variation we find O = 0.334 + 0.119z. As we can see, the profile variations tilt the structure on the outboard 
mid-plane. This finding agrees qualitatively with other published results for full global simulations of linear 
ITG modes [29]. Our interpretation is that this is a consequence of profile variations in the equilibrium. 


3.3 Calculations with Profile Variations and Sheared Toroidal Flow 

We now consider the influence of a small toroidal flow shear on stability and the reconstructed mode 
structures. This has been achieved by moving into the frame rotating toroidally with the plasma surface at 
x = 0, and by introducing the following Doppler shift into 2 : 


D 0 O,p) -+ £lo(x,p) - ntyx = 


vo(x,p) - ntt^x 


+ Z7o0,p), 


( 8 ) 


where is constant and sets the flow shearing rate. Thus fl'^x represents the toroidal rotation frequency 
(normalised to c s /a) of the magnetic flux surfaces measured relative to the mode rational surface at x = 0 
where the reconstructed global mode sits. We have considered —0.06 < < +0.06. Note that we assume 

a purely toroidal flow here, as poloidal flows are expected to be damped in a tokamak plasma. From Eq. ( 8 ) 
it is clear that the flow shear influences the coefficients of Eq. (5). Flow shear will therefore influence 
global modes in a similar way to the variations in other profiles, which we will now quantify. 

Depending on the sign of the flow shear the reconstructed global mode shifts upward or downward 
relative to the outboard mid-plane with a structure in the poloidal plane that is very similar to that in 
Figure 7. For a mode peaking at 6 = 0 in the absence of flow shear, increasing flow shear tilts the mode 
structure on the outboard mid-plane and lowers the linear growth rate. Figure 8 shows how the growth rate 
varies with flow shear for the different profile variations considered earlier in Figure 4 and Figure 7. Keeping 
only the symmetric Lt and L n profiles the growth rate curve is symmetric about = 0, shown by curve 
(e) of Figure 8 . However, when other profile variations are taken into account, an asymmetry is introduced 
into the dependence of growth rate on the sign of as shown in Figure 8 (for curves a-d). Only in the 
case where we have purely symmetric profiles about x = 0 is the growth rate maximised at = 0. In 
the other cases the peak in growth rate occurs at non-zero flow shear. We can explain this as follows. We 
expect this maximally unstable isolated mode to exist when both u)${x,p) and 70 (x,p) are stationary at the 
same position. From symmetry we anticipate duJo/dp\ p =o = $7o/dp|p=o = 0, so let us consider p = 0. Then 
Doppler-shifting our expression for Qq( x,p) given in Eq. (5) and differentiating with respect to x, we find 


duo 

dx 


= a, J 1 ) - nQ'j. + 2 a! 2 ); 


and 


p =0 

d7Q 

dx 


= + + 2a®* 


p =0 


(9) 


( 10 ) 


where affi = a k-‘ ^ an d th e subscripts r and i indicate the real and imaginary component respectively. We 

require d^fo/dx = 0 , which then provides an equation for the mode’s radial position, x = xo = —a[^/2a[ 2 \ 
Substituting this into Eq. (9) we find the critical shearing rate, Q^ m for which duio/dx = 0 at this same 
value of x = xq:~ 


O' — — 
n 


„( 2 )Jl) 

71) a r a i 
( 2 ) 


a 


(ii) 


2 While this method can in principle handle an arbitrary experimental profile of toroidal rotation, for the simple example 
here we assume that the toroidal flow varies linearly in x. 
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Thus, we expect the maximally unstable isolated mode to exist for this critical shearing rate, = Q^ m , 
but centred on x = xq rather than x = 0. Considering the full profile variation case (curve-d), our numerical 
results show the growth rate peaks at « —0.017 with a value of 7 = 0.136, which is ~ 10% higher than 
the growth rate for zero flow shear (7 = 0.119). For this case = —0.668+ 0.269z and a= 0.361 — 6 . 68 z. 
Substituting these values into Eq. (11) provides Q^ m = —0.0168, which is in good agreement with curve-d 
of Figure 8 . 

We can compare our flow shear results with the global simulations performed in Ref [29]. This is 
complicated as we employ a toroidal flow, while the simulations of Ref [29] employ an E x B flow, which 
is almost poloidal. Nevertheless, if parallel flows have a negligible impact, the two can be related by a 
geometric factor. We can factor out this geometric factor by considering the ratio of the flow shear which 
maximises the growth rate to the value required to stabilise the mode. Our result of —0.38 then agrees very 
well with that of Ref [29], which is —0.39. 

Figure 9 presents the scaling of fiW and offset in the poloidal angle po with toroidal mode number n. As 
expected from Eq. (11) scales with 1/n (or p*). Furthermore, we find that at large n \po\ = a + bn ~ a , 

where for our profile choice a = 0.041, b = 1.88 and a = 0.925. We have also investigated the effect of 
asymmetry in the growth rate spectrum on the reconstructed global mode structure, and this is illustrated 
in Figure 10. For = 0 the structure is already tilted, but increasing flow shear in the negative direction 
acts to re-align the mode radially and for a critical value of flow shear, « —0.017, the effect of the 

profile variation is completely compensated, allowing an isolated mode again to form with largest growth 
rate, 7 =Max[ 7 o(x,p)]. Note the mode is radially shifted slightly relative to x = 0 at this value of This 
is consistent with the above analysis. Increasing flow shear even further, beyond the critical value, tilts the 
mode structure in the opposite direction and lowers its linear growth rate again. These results, obtained 
purely from solutions of GS 2 and the higher order theory, are again in good qualitative agreement with 
global calculations of linear electrostatic ITG modes presented in Ref [29] . 

To summarise, for realistic and experimentally relevant cases where we take the profile variations into 
account, we do not in general expect to find pure isolated modes. Isolated modes can form only in special 
radial locations where the equilibrium profiles produce a stationary point in f2o(a?,p). However, making 
adjustments to one equilibrium profile while the others are fixed, can produce the required stationary point 
and lead to the onset of the isolated mode, as arises in the above example for a critical toroidal flow shear 
equal to 0 ^ m . 

4 Conclusion 

In this work we have reconstructed the 2D global mode structure and the global growth rate for linear 
electrostatic ITG modes, using local solutions from a gyrokinetic code (GS 2 ) and higher order ballooning 
theory. This approach, which is solid provided that equilibrium quantities vary slowly across rational 
surfaces, has provided additional insight into the physics of global simulations of linear microinstabilities 
in tokamak plasmas. Our first investigations used radial profiles for the mode drives that were peaked 
and symmetric about x = 0 , and we held all other equilibrium profiles constant; this results in the local 
complex mode frequency, f2o(aJ,p), having a stationary point at x = 0. This condition produces a special 
class of mode, known as the “isolated mode”, that peaks at the outboard mid-plane with a large growth 
rate, 7 ~Max[ 7 o(#,p)]. These results are in very good qualitative agreement with the simplified fluid model 
of ITG modes presented in Ref [19]. Introducing radial variation into other equilibrium profiles, we show 
that the radial position of the stationary points in both local frequency, ujq , and growth rate, 70 , become 
shifted with respect to each other. In this case, the reconstructed global mode becomes less unstable and 
shifts poloidally away from the outboard mid-plane. 

Toroidal flow shear, introduced as a Doppler shift in the real frequency, also influences the global mode. 
Starting from the conditions of an isolated mode, with drive profiles peaked and symmetric about x = 0 
and with no other profile variations, adding a constant flow shear is always found to be stabilising. When 
other profile variations are included, flow shear can be destabilising when the flow shear counteracts the 
tilting of the mode structure at the outboard mid-plane that is induced by the other profile variations. This 
results in an asymmetry in the growth rate as a function of flow shear about = 0 , which is in qualitative 
agreement with previous global gyrokinetic calculations [29] . Moreover, flow shear is also found to shift the 
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mode radially. For a critical flow shear (or a critical toroidal mode number for a given flow shear - see 
Eq. (11) ) the isolated mode can exist even with arbitrary profiles. 

In this paper we have focussed mainly on studying special surfaces where fast growing isolated modes 
arise because of a peaked drive that results in a stationary point in Hq (x,p). For an arbitrary equilibrium 
the drive for a particular mode number may indeed have local maxima on some special surfaces, which can 
be located using local gyrokinetic solutions to find stationary points in Qq(x,p). It is this strongly growing 
isolated mode that will typically be captured by global, initial value gyrokinetic codes. Nevertheless, away 
from these special surfaces we do not expect to find high growth rate isolated modes, but slower growing 
general modes will predominate. (We note that local codes cannot distinguish between isolated and general 
modes without including higher order corrections.) This might have important implications for quasilinear 
transport models, and for flow generation mechanisms, though other effects associated with nonlinearities 
and turbulence will also be important: fully nonlinear global simulations are needed to assess the relative 
importance of the linear physics that is described here. In this paper we have exploited the higher order 
analysis to obtain the global mode structure and frequency of the stronger isolated mode. Future work will 
assess the impact of the slower growing general modes on transport in other regions of the plasma away 
from stationary points in Ho 

Finally, we point out that the procedure used in this work is quite general and can be used to explore 
more realistic tokamak equilibria and more complicated instabilities, including the effect of shaping (e.g. 
elongation and triangularity) and electromagnetic modes with kinetic electrons. Future work will extend 
our study to explore these effects. 
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